home *** CD-ROM | disk | FTP | other *** search
/ Monster Media 1996 #15 / Monster Media Number 15 (Monster Media)(July 1996).ISO / prog_c / cuj0696.zip / DWYER.ZIP / SPECTRAL.TST / SPECTRAL.CX < prev    next >
Text File  |  1996-03-23  |  8KB  |  318 lines
  1. /* ============ */
  2. /* spectral.cx    */
  3. /* ============ */
  4. #include <math.h>
  5. #include <stdio.h>
  6. #include <stdlib.h>
  7. #include <string.h>
  8. #include <assert.h>
  9. #include <xtendefs.h>
  10. #include <specdefs.h>
  11.  
  12. #include <calcmert.cx>
  13. #include <euclstep.cx>
  14. #include <getctlda.cx>
  15. #include <getgcd.cx>
  16. #include <prtspdat.cx>
  17. #include <putxflt.cx>
  18. #include <putxfrac.cx>
  19. #include <putxint.cx>
  20. #include <xdot.cx>
  21. #include <xmatprt.cx>
  22. /* ==================================================================== */
  23. /* MAIN - Performs Spectral Test on multiplier of CRNG            */
  24. /* ==================================================================== */
  25. void
  26. main()
  27. {
  28.     short   HowMany, InputOK, MaxDimen;
  29.     USHORT  Modulus[NE], Mult[NE * 2];
  30.     USHORT  Tmp1[NE], Tmp2[NE];
  31.  
  32.     InputOK = GetControlData(Mult, Modulus, &MaxDimen, &HowMany);
  33.  
  34.     if (InputOK)
  35.     {
  36.     USHORT    XH[NE], XHP[NE], XP[NE], XPP[NE],
  37.         XR[NE], XS[NE],  XU[NE], XV[NE];
  38.     USHORT    MSqd[NE];        /* Modulus Squared      */
  39.     USHORT    XMerit[NE];        /* Holds Accuracy Results */
  40.  
  41.     printf("\n%22d%41s", HowMany, "Number of Multipliers");
  42.     printf("\n%20s", "");
  43.     PutXInt(Mult, 20);
  44.     printf("  Multiplier");
  45.     if (HowMany == 1)
  46.     {
  47.         /* Do Nothing */
  48.     }
  49.     else
  50.     {
  51.         /* ------------------------------- */
  52.         /* Print Label on First Multiplier */
  53.         /* ------------------------------- */
  54.         printf(" # 1\n%20s", "");
  55.  
  56.         /* ----------------------- */
  57.         /* Print Second Multiplier */
  58.         /* ----------------------- */
  59.         PutXInt(Mult + NE, 20);
  60.  
  61.         /* ----------- */
  62.         /* Print Label */
  63.         /* ----------- */
  64.         printf("  Multiplier # 2");
  65.     }
  66.  
  67.     printf("\n%20s", "");
  68.     PutXInt(Modulus, 20);
  69.     printf("  Modulus\n%21s", "");
  70.     printf("%-19d  Maximum Dimension of Test\n", MaxDimen);
  71.  
  72.     if (HowMany == 1)
  73.     {
  74.         /* -------------- */
  75.         /* S1: Initialize */
  76.         /* -------------- */
  77.  
  78.         XCOPY(Mult, XH);        /* h <- a    */
  79.         XCOPY(Modulus, XHP);    /* m <- h'      */
  80.         XCOPY(EONE, XP);        /* p <- 1    */
  81.         XCOPY(EZERO, XPP);        /* p' <- 0      */
  82.         XCOPY(Mult, XR);        /* r <- a    */
  83.         XSQR(Mult, XS);
  84.         XADD(XS, EONE, XS);        /* S <- 1 + a^2 */
  85.  
  86.         /* ------------------ */
  87.         /* S2: Euclidean Step */
  88.         /* ------------------ */
  89.         EuclStep(XH, XHP, XP, XPP, XS, XU, XV);
  90.  
  91.         /* -------------- */
  92.         /* S3: Compute V2 */
  93.         /* -------------- */
  94.         XSUB(XU, XH, XU);        /* u -= h    */
  95.         XSUB(XV, XP, XV);        /* v -= p    */
  96.         XSQR(XU, Tmp1);        /* u^2        */
  97.         XSQR(XV, Tmp2);        /* v^2        */
  98.         XADD(Tmp1, Tmp2, Tmp2);    /* u^2 + v^2    */
  99.         if (XGT(XS, Tmp2))        /* if(s > Tmp2) */
  100.         {
  101.         XCOPY(Tmp2, XS);    /* s <- u^2+v^2 */
  102.         XCOPY(XU, XHP);        /* h' <- u      */
  103.         XCOPY(XV, XPP);        /* p' <- v      */
  104.         }
  105.     }
  106.     else
  107.     {
  108.         XSQR(Modulus, XS);
  109.         XCOPY(XS, MSqd);
  110.     }
  111.  
  112.     if (HowMany == 1)
  113.     {
  114.         CalcMerit(XS, 2, Modulus, XMerit);
  115.     }
  116.     else
  117.     {
  118.         CalcMerit(XS, 2, MSqd, XMerit);
  119.     }
  120.  
  121.     PrintSprectalData(XS, XMerit, 2);
  122.  
  123.     if (MaxDimen > 2)
  124.     {
  125.         short   I, J, K, NxtDimen;
  126.         USHORT *Umat, *Vmat;
  127.         USHORT  R11[NE], R12[NE], R21[NE], R22[NE];
  128.         USHORT  Vij[NE], Vjj[NE];
  129.         USHORT  XA[NE], XB[NE];
  130.         USHORT *X, *Y, *Z;
  131.  
  132.         /* ----------------- */
  133.         /* Initialize Arrays */
  134.         /* ----------------- */
  135.         INIT_VEC(X, MaxDimen);
  136.         INIT_VEC(Y, MaxDimen);
  137.         INIT_VEC(Z, MaxDimen);
  138.  
  139.         INIT_MAT(Umat, MaxDimen);
  140.         INIT_MAT(Vmat, MaxDimen);
  141.  
  142.         /* -------------------------------- */
  143.         /* Initialize U and V Matrices Etc. */
  144.         /* -------------------------------- */
  145.         if (HowMany == 1)
  146.         {
  147.         INIT_UMAT(Umat, MaxDimen);
  148.         INIT_VMAT(Vmat, MaxDimen);
  149.         }
  150.         else
  151.         {
  152.         XCOPY(Mult, XA);
  153.         XCOPY(Mult + NE, XB);
  154.         INIT_MAT_2(Umat, Modulus, MaxDimen);
  155.         INIT_MAT_2(Vmat, EONE, MaxDimen);
  156.         INIT_R_MAT(R11, R12, R21, R22);
  157.         }
  158.  
  159.         MATPRT(Umat, 2, MaxDimen, "Initial U Matrix");
  160.         MATPRT(Vmat, 2, MaxDimen, "Initial V Matrix");
  161.  
  162.         /* --------------------- */
  163.         /* S4: Advance Dimension */
  164.         /* --------------------- */
  165.         for (NxtDimen = 3; NxtDimen <= MaxDimen; ++NxtDimen)
  166.         {
  167.         USHORT    Utt[NE];
  168.  
  169.         INIT_VNEW(Vmat, NxtDimen, Modulus, MaxDimen);
  170.  
  171.         if (HowMany == 1)
  172.         {
  173.             XMULT(XR, Mult, Tmp1);    /* r * a */
  174.             XMOD(Tmp1, Modulus, XR);    /* r <- (ra)mod m */
  175.  
  176.             INIT_UNEW(Umat, NxtDimen, XR, MaxDimen);
  177.  
  178.             FINISH_UV(Umat, Vmat, NxtDimen, Modulus, MaxDimen, XR);
  179.         }
  180.         else
  181.         {
  182.             SET_NEXT_RMAT(R11, R12, R21, R22, XA, XB, Modulus);
  183.             SET_NEXT_UROW(Umat, NxtDimen, R12, R22, MaxDimen);
  184.  
  185.             COMPLETE_UV(Umat, Vmat, NxtDimen, Modulus, MaxDimen,
  186.             R12, R22);
  187.         }
  188.         MATPRT(Umat, NxtDimen, MaxDimen, "U Matrix At Step S5");
  189.         MATPRT(Vmat, NxtDimen, MaxDimen, "V Matrix At Step S5");
  190.  
  191.         DOT(NxtDimen, Umat, NxtDimen, Umat, NxtDimen, MaxDimen, Utt);
  192.  
  193.         XMIN(Utt, XS, XS);
  194.         J = 1, K = NxtDimen;
  195.  
  196.         do
  197.         {
  198.             /* ------------- */
  199.             /* S5: Transform */
  200.             /* ------------- */
  201.             DOT(NxtDimen, Vmat, J, Vmat, J, MaxDimen, Vjj);
  202.  
  203.             for (I = 1; I <= NxtDimen; ++I)
  204.             {
  205.             if (I != J)
  206.             {
  207.                 DOT(NxtDimen, Vmat, I, Vmat, J, MaxDimen, Vij);
  208.                 XCOPY(Vij, Tmp1);
  209.                 XABS(Tmp1);
  210.                 XRMULT(Tmp1, ETWO);
  211.                 if (XGT(Tmp1, Vjj))
  212.                 {
  213.                 TRANSFORM(Vmat, I, Umat, J, Vij, Vjj,
  214.                     NxtDimen, MaxDimen);
  215.                 K = J;
  216.                 }
  217.             }
  218.             }
  219.  
  220.             /* --------------------- */
  221.             /* S6: Examine New Bound */
  222.             /* --------------------- */
  223.             if (K == J)
  224.             {
  225.             USHORT    Ujj[NE];
  226.  
  227.             DOT(NxtDimen, Umat, J, Umat, J, MaxDimen, Ujj);
  228.             XMIN(XS, Ujj, XS);
  229.             }
  230.             /* ------------- */
  231.             /* S7: Advance J */
  232.             /* ------------- */
  233.             if (J == NxtDimen)
  234.             {
  235.             J = 1;
  236.             }
  237.             else
  238.             {
  239.             ++J;
  240.             }
  241.         }
  242.         while (J != K);
  243.  
  244.         MATPRT(Umat, NxtDimen, MaxDimen, "U Matrix At End of S7");
  245.         MATPRT(Vmat, NxtDimen, MaxDimen, "V Matrix At End of S7");
  246.  
  247.         /* ----------------------- */
  248.         /* S8:    Prepare for Search */
  249.         /* ----------------------- */
  250.         CLR_VECTOR(NxtDimen, X);
  251.         CLR_VECTOR(NxtDimen, Y);
  252.  
  253.         for (J = 1; J <= NxtDimen; ++J)
  254.         {
  255.             XDIV(XS, Modulus, Tmp1);
  256.  
  257.             DOT(NxtDimen, Vmat, J, Vmat, J, MaxDimen, Z + XVL(J));
  258.             XRDIV(Z + XVL(J), Modulus);
  259.             XRMULT(Z + XVL(J), Tmp1);
  260.             XFLOOR(Z + XVL(J), Z + XVL(J));
  261.             XSQRT(Z + XVL(J), Z + XVL(J));
  262.             XFLOOR(Z + XVL(J), Z + XVL(J));
  263.         }
  264.  
  265.         VECPRT(Z, NxtDimen, "Z Vector at Step S8");
  266.  
  267.         /* ------------------------------------------------ */
  268.         /* S9 - S11 - The for control statement is Step S11 */
  269.         /* ------------------------------------------------ */
  270.         for (K = NxtDimen; K >= 1; --K)
  271.         {
  272.             /* -------------- */
  273.             /* S9: Advance Xk */
  274.             /* -------------- */
  275.             if (!XEQUAL(X + XVL(K), Z + XVL(K)))
  276.             {
  277.             XRADD(X + XVL(K), EONE);
  278.             INCREASE_Y(Y, Umat, K, NxtDimen, MaxDimen);
  279.             VECPRT(Y, NxtDimen, "Increased Y: ");
  280.             do
  281.             {
  282.                 /* -------------- */
  283.                 /* S10: Advance K */
  284.                 /* -------------- */
  285.                 ++K;
  286.                 if (K > NxtDimen)
  287.                 {
  288.                 XdotProd(NxtDimen, Y, Y, Tmp1);
  289.                 XMIN(XS, Tmp1, XS);
  290.                 }
  291.                 else
  292.                 {
  293.                 XCOPY(Z + XVL(K), X + XVL(K));
  294.                 XNEG(X + XVL(K));
  295.                 DECREASE_Y(Y, Umat, K, Z + XVL(K),
  296.                     NxtDimen, MaxDimen);
  297.                 VECPRT(Y, NxtDimen, "Decreased Y: ");
  298.                 }
  299.             }
  300.             while (K <= NxtDimen);    /* end of Step S10   */
  301.             }            /* end of if  - Step S9      */
  302.         }            /* end of for - Steps S9-S11 */
  303.  
  304.         if (HowMany == 1)
  305.         {
  306.             CalcMerit(XS, NxtDimen, Modulus, XMerit);
  307.         }
  308.         else
  309.         {
  310.             CalcMerit(XS, NxtDimen, MSqd, XMerit);
  311.         }
  312.  
  313.         PrintSprectalData(XS, XMerit, NxtDimen);
  314.         }                /* end of for - Step S4 */
  315.     }
  316.     }
  317. }
  318.